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Abstract 

We study the effect of global error control in the numerical solution 
of Hamiltonian systems. In particular, we apply the RKQ algorithm 
in the numerical solution of a Hamiltonian system. This algorithm 
is designed to provide stepwise control of both local and global er- 
ror. A test problem demonstrates the error control features of RKQ. 
Good results are obtained, despite the fact that explicit Runge-Kutta 
methods have been used in RKQ, rather than symplectic Runge-Kutta 
methods. This simply emphasizes the value of stepwise global error 
control, as per the RKQ algorithm. 

1 Introduction 

Global error control in numerical solutions of initial-value problems is of 
paramount importance. In this paper, we consider the effect of global error 
control on the numerical solution of Hamiltonian systems. We solve a test 
problem using the RKQ algorithm, which is designed to achieve stepwise 
global error control. Significantly, we use explicit Runge-Kutta methods in 
RKQ, as opposed to implicit symplectic Runge-Kutta methods, which are 
often preferred for Hamiltonian systems. 
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2 Relevant Concepts, Terminology and No- 
tation 



2.1 Hamiltonian Systems 

If the total energy of a physical system is given in the form of the Hamiltonian 

H (p, q) , 

where the canonical coordinates p and q are the momentum and position, 
respectively, then the evolution of the system is given by 

dq dH dp dH 
dt dp ' dt dq 

The Hamiltonian H (p, q) is a first integral of flTJ, meaning that it is constant 
for all t. For an appropriate set of initial values, ([T]) is an initial-value problem 
(IVP). Since Hamiltonian systems are nonlinear, their solution is usually 
obtained numerically - using, for example, a Runge-Kutta (RK) method. 

For ease of presentation we consider here a two-dimensional Hamiltonian, 
rather than the more general case 

H(pu ...,p n ,qi,---,qn), 

but this restriction will not affect our discussion. 



2.2 Runge-Kutta Methods 

Runge-Kutta methods are very well-known, and the reader is referred to the 
extensive literature. For our purposes, an RK method applied to ([T]) has the 
form 

*+ J ] = [*]+ h i+1 F (tt, \ ® 1 , h i+1 ) . (2) 

. Pi+i J L Pi \ v L ^ J J 

In (j2j), i denotes discrete points along the t-axis, so that q~i and are numer- 
ical solutions at tf, the stepsize ftj+i = ti+\ — tf, and F is a vector function 
associated to the RK method under consideration. Initial values of p and q 
are specified at to- 

We define the local and global errors of an RK method as 



local error: e 



i+l 



Qi 
Pi 



+ h i+1 F t 



q, 
pi 



! hi+1 



q-i+i 
Pi+i 



2 



and 



global error: A i+ i = 




Qi+i 
Pi+i 



Qi+i 
Pi+i 



(3) 



Note the use of the exact values of q and p in the definition of the local error. 
If the RK method is of order r (denoted RKr), we have 



where h is representative of the stepsizes (for example, h could be taken as 
the average stepsize along the discretized t-axis). 

Usually, symplectic RK methods are used to solve Hamiltonian systems 
PP. Such methods have the property that the Hamiltonian arising from the 
numerical solution is bounded in the sense that it does not drift away from 
the exact value; rather, it exhibits small oscillations in the vicinity of the 
exact value. Also, symplectic RK methods reproduce closed orbits in the 
(q,p) phase space, as expected when periodic solutions are present. The 
symplectic property of Hamiltonian systems refers to the invariance of the 
differential 2-form dp A dq, a feature that is respected by symplectic RK 
methods (hence their name). Practically speaking, however, it is the first 
two properties listed here - essentially constant numerical Hamiltonian and 
closed phase space trajectories - that are of primary interest. 

The disadvantage of symplectic RK methods is that they are implicit (a 
particular characteristic of the function F), and this requires the solution of 
a nonlinear system of equations at each node ti, which can be computation- 
ally expensive. Explicit RK methods do not require the solution of such a 
nonlinear system, but neither are they symplectic. 



We will not discuss KKrvQz in detail here; the reader is referred to our 
previous work where the algorithm has been discussed extensively [2j,|2]- It 
is sufficient to state that KKrvQz uses RKr and KKv to control the local 
error via so-called local extrapolation, while simultaneously using KKz to 
keep track of the global error in the RKr solution (we have z 3> r, v). Such 
global error arises from the propagation of the KKv global error in the RKr 
method, as a consequence of local extrapolation. KKrvQz is designed to 
estimate the various components of the global error in RKr and KKv at 
each node ti and, when the global error exceeds a user-defined tolerance, a 
quenching procedure is carried out. This simply involves replacing the RKr 
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2.3 The KKrvQz Algorithm 
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and RKi> solutions with the much more accurate KKz solution, whenever 
necessary, so that the RKr and KKv global errors do not accumulate beyond 
the desired tolerance. 



3 The Effect of Global Error Control 

From (l3l) we have 



Pi+i 



Qi+l + $q,i+l 
Pi+1 + Op,i+l> 



This gives 
H(p i+1 ,q i+l ) 



rji \ . dH (p i+ i,g i+ i) 
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dH (p i+1 ,q i+1 ) 
dq 



where we have ignored higher-order terms. Now, say 5 is an upper bound on 
the magnitude of the global errors. Hence, 



\H (p i+1 ,q i+1 ) - H (p i+ll q i+ i)\ < 
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We see that the bound on the error in the numerical Hamiltonian is propor- 
tional to the bound S. 

For a point on the trajectory in phase space, we have 



\\(q i+ i,p i+1 ) - (q i+1 ,p i+1 ) 



where ||---|| is any norm suitable for determining distances in the phase 
space. Hence, the bound on the trajectory error is proportional to 5. For the 
Euclidean norm and the two-dimensional case considered here, we have 

\\(q i+1 ,p i+1 ) - (q i+h p i+ i)\\ ^ y/25. 

The implication of the above analysis is obvious: if we apply KKrvQz 
to the problem, with a suitable tolerance of 5, then we would generate solu- 
tions p and q for which the error bounds on the numerical Hamiltonian and 
the phase-space trajectories are acceptably small. Furthermore, this can be 
achieved using explicit RK methods in the KKrvQz algorithm, as opposed 
to the more computationally intensive symplectic RK methods. 
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4 Numerical Example 



As an example, we consider the Hamiltonian 



H 



( 




cosg 



which yields 



dq 

~di 

dp 

~dt 
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This is the same example considered in Hairer et al [T] . We use initial values 



and we integrate over t G [0, 4000] . 

We solve the system using RK34Q8 (r = 3, v — 4, z — 8) , with a tolerance 
of S = 10 -6 on both the local and global error. For comparison, we solve 
the system via local extrapolation only with RK3 and RK4, also subject to 
a tolerance (on the local error) of 5 = 10 -6 . Results are shown in the figures 
following the references. The various explicit RK methods used here are the 
same as those referenced in [2], [3]. We use the notation RK34 to indicate the 
local extrapolation algorithm using RK3 and RK4. 

In Figure 1, the top two plots show the first few periods of the solution 
numerical p(t) and q(t) , demonstrating their periodic character. Consistent 
with this periodicity is the closed trajectory in phase space, shown in the third 
plot in Figure 1. The numerical Hamiltonian is shown in the fourth plot, for 
RK34Q8 and RK34. Clearly, the latter exhibits a generally monotonic drift 
from the exact value of 0.8, while the former is essentially constant with slight 
oscillations. For the sake of clarity we have not shown all the data points in 
this plot; there are some 96000 nodes on [0,4000] in the computation, and 
we show relatively few - sufficient, nonetheless, to exhibit the salient features 
of the calculation (this also holds for subsequent plots in Figures 2-4). The 
drift in H is slight - only about 3 x 10~ 6 over the entire interval of integration 
- but it is definite. Given that H should be invariant, however, the result 
obtained with RK34Q8 should be preferred. 

Figure 2 shows global errors in p (t) and q (t) , determined from the dif- 
ference of the RK34 solution and the RK8 solution. Clearly, the errors in 
the RK34 algorithm grow as the integration proceeds, achieving maximal 
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values of ~ 0.02. On the other hand, the RK34Q8 algorithm, thanks to the 
high-order quenching device, gives a solution that is bounded by 5 = 10~ 6 , as 
desired. This bound is indicated by the horizontal line labelled 'tolerance'. 
The maximal errors in this case are 9.998 x 10~ 7 in q (t) , and 8.54 x 10~ 7 in 
P(t) ■ 

In Figure 3 we show the trajectory error, determined using the Euclidean 
norm, for both algorithms. Clearly, the trajectory error grows for RK34, 
consistent with growth of global error in the solutions seen in Figure 2. For 
RK34Q8, the trajectory error is bounded by a/2 x 10~ 6 , as expected. This 
bound is indicated by the horizontal line labelled 'upper bound'. 

For completeness, we show the estimated global error in the RK8 solution 
in Figure 4 (errors for both p and q are shown). This has been computed 
with (see 0,13, @]) 

= £ i+ \ + (I 2 + h i+ iF q ^ p ) Aj 

~ + (I2 + hi+lfq,p) Aj. 



Here, F 9iP is the Jacobian of the function F in 
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I > is the 2x2 identity matrix, and we assume 



), f ?iP is the Jacobian of 








Also, we have used Richardson extrapolation to estimate the local error £j + i 
in the RK8 solution [2], [5]. For a reasonably small stepsize (~ 0.05 in this 
problem) and a high-order RK method, we expect this to be a good estimator. 
The largest error in either component of the RK8 solution was 1.86 x 10~ 12 . 
This means that any global error in the RK8 solution would not 'contaminate' 
the RK34Q8 error control procedure, since the tolerance of 5 = 10 -6 is 
considerably larger than 1.86 x 10 -12 . We note here that it is feasible to use 
a method of higher order than RK8 to estimate the global error in the RK8 
solution, even though this would probably increase the computational effort. 
At the time of writing, however, we did not have access to such a method. 
Nevertheless, in this regard one might consider transforming the problem 
into a second-order problem, making it suitable for Nystrom integration, as 
we considered in [6]. This would enable the use of the (10,12) embedded 
Nystrom pair due to Dormand et al [7]. 
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4.1 Relative Error Control 

In this paper, we have only considered absolute error control, since the mag- 
nitude of the solution was ~ 1. If the magnitude of the solution had been 
much larger than unity, it would have been better to implement relative er- 
ror control. We will not discuss this in detail; a thorough account has been 
given in [3j, wherein we generalized the RKQ algorithm. It suffices to say 
that at each node, for relative error control in the problem considered here, 
the tolerance would have the form 

8 = mm {max {5 A , 5 R \qi\} , max {5 A , 5 R \pi\}} , 

where 5a and 8r are user-defined (the presence of 8a caters for those situa- 
tions where qt and/or are very close to zero). 

5 Conclusion 

We have explored global error control in the numerical solution of Hamilto- 
nian systems. A theoretical analysis shows that if the error in the solution 
is bounded, then errors in quantities such as the numerical Hamiltonian and 
phase space trajectories are also bounded. We have considered the use of the 
RKQ algorithm, with explicit Runge-Kutta methods, as our choice of nu- 
merical integrator, since RKQ is specifically designed to achieve global error 
control. A test problem has demonstrated the expected results. In addition 
to the bounding of the error in the solution, we also observe an essentially con- 
stant numerical Hamiltonian and a bounded error in the numerical trajectory 
in phase space. This contrasts Hamiltonian drift and unbounded trajectory 
deviation normally associated with the use of explicit Runge-Kutta methods 
in solving Hamiltonian systems. It is our contention that, even though RKQ 
utilizes explicit Runge-Kutta methods, stepwise global error control via RKQ 
leads to results with a similar quality to those that would be obtained using 
symplectic methods. 
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